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O ■ Abstract 
m 

We present data analysis methods used in detection and the estimation of parameters of gravita- 

Qs, ' tional wave signals from the white dwarf binaries in the mock LISA data challenge. Our main focus 

O ■ 

is on the analysis of challenge 3.1, where the gravitational wave signals from more than 6 x 10 7 
^ . Galactic binaries were added to the simulated Gaussian instrumental noise. Majority of the signals 

at low frequencies are not resolved individually. The confusion between the signals is strongly 
reduced at frequencies above 5 mHz. Our basic data analysis procedure is the maximum likelihood 
detection method. We filter the data through the template bank at the first step of the search, 
then we refine parameters using the Nelder-Mead algorithm, we remove the strongest signal found 
and we repeat the procedure. We detect reliably and estimate parameters accurately of more than 
ten thousand signals from white dwarf binaries. 

PACS numbers: 95.55.Ym, 04.80.Nn, 95.75.Pq, 97.60.Gb 
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I. INTRODUCTION 



The Galaxy contains a numerous population of the ultra-compact binaries which have 
orbital periods shorter than one hour. The observations of those binaries carry important 



astrophysical in 
their evolution 



brmation about the internal stellar structure, formation of binaries and 
1]. Due to the short period and proximity, the ultra-compact Galactic 
binaries will be an important source of gravitational waves (GW) for the future space borne 
interferometer LISA. LISA is planned to be launched in the next decade jointly by ESA and 
NASA. The bandwidth of the LISA detector is expected to be from 0.1 mHz to 100 mHz. In 
this frequency range we expect around 6 x 10 7 ultra-compact binaries. These binaries will 
be dominated by the population of the white dwarf binaries. The number of the observed 
ultra-compact binaries will be so large below 3 mHz that they are not individually resolvable 
and form a cyclo-stationary background which dominates over the instrumental noise above 
0.1 mHz {2I, yj|. The number of binaries drops significantly above 7-8 mHz. 

The most common sources are white-dwarf /white-dwarf binaries emitting gravitational 
wave signals of nearly constant frequency and amplitude. We also expect to observe few 
white-dwarf /neutron star binaries. The binaries could be of two major types: 

(i) Detached, separated white-dwarf /white-dwarf binaries whose evolution is driven by 
radiation reaction. They are the end points of many binary evolution scenarios. The gravi- 
tational wave carry information about the mass of the binary and the distance. 

(ii) Interacting binaries. Those are close systems with a significant tidal interaction 
and/or with the Roche lobe overflow. In those systems the gravitational radiation reaction 
competes against mass transfer and the orbital period can either increase or decrease. Cur- 
rently there are 22 known accreting binaries, so called AM CVn stars, with periods between 
5.4 and 65 min |4J. The radiation from those binaries fall into the LISA band and they will 
serve as verification binaries 15 1 . 

n 

It was demonstrated [6[ that one can detect and remove a few tens of thousand of those 
signals. The resolved systems will provide the map of the compact binaries in the Galaxy 
and will allow us to constrain the evolutionary pathways of those systems. 

A series of mock LISA data challenges (MLDC) was organized in order to foster develop- 
ment of LISA data analysis algorithms and to compare the performance of different methods 
[7 10 1 . The simulated Galaxy consists of only white dwarf binaries (of both types detached 
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and interacting) in circular orbit and it is based on the population synthesis described in 
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12] . It is expected that white dwarf binaries will largely dominate over more heavy bina- 



3e gene rated 
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14j. The 



ries like neutron stars and/or black holes. Eccentric white dwarf binaries could 
in globular clusters but their number is a tiny fraction of the total population 
gravitational wave signals produced by white dwarf Galactic binaries span the whole LISA 
band, they stand above the instrumental noise starting at 0.1 mHz and propagate all the 
way up to few tens of mHz. 

Let us give a brief overview of currently available methods. The fully coherent methods 
employing matched filtering techniques can be split in two groups: stochastic search and 
grid-based search. 

The stochastic search does not map uniformly the whole parameter space, instead it 
concentrates on the regions with high likelihood. The first type of the stochastic search 
is suggested in [ill and it is based on the genetic algorithm. The genetic algorithm is 
an optimization method which evolves the set of templates (a colony of organisms) in the 
direction of increasing likelihood (improving fitness of organisms) using a certain rules (selec- 
tion, breeding, mutation). Another stochastic method is based on constructing the (Markov) 
chains using Metropolis-Hastings acceptance/rejection rule. The Bayesian methods are pow- 
erful tools to get posterior probability distribution function. A pure MCMC (Markov chain 
Monte-Carlo) algorithm does not perform well due to presence of quite strong (and well sep- 
arated in the parameter space) secondary maxima in the likelihood. Two currently available 
Bayesian algorithms differ mainly in the way they explore the parameter space. BAM algo- 
rithm suggested in 16] uses multiple proposal distributions reflecting possible correlations 
in the parameter space in combination with the simulated annealing. At the search stage 
(the search for the global maximum in the parameter space) the Markovian properties of 
the chain are quite often not respected, and, once the search is completed, the sampling 
stage starts during which the classical MCMC algorithm is used. The key feature of BAM is 
blocking: the search is conducted in the several frequency bands each split in several small 
blocks. The algorithm steps through these blocks updating all sources within a given band 
simultaneously. After all blocks have been updated, they are shifted by one-half the width 
of a blocks for the next round of updates to eliminate boundary effect. The results from 
the different frequency bands are glued together using the overlapped buffer zones. Another 



MCMC-based method is described in 



17] . The authors suggest to use delayed rejection 
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MCMC to explore the parameter space, find the global maximum and sample the posterior 
distribution. The basic idea behind "delayed rejection" is an extended acceptance/rejection 
rule: we allow several jumps before rejecting/accepting the next point in the chain, and each 
following trial jump learns the property of the parameter space explored by the previous 
jumps. An additional advantage of the "delayed rejection" is that the variance of an esti- 
mate made from a chain using delayed rejection is always smaller than that produced with 
a standard Markov chain. 

The grid-based search maps the whole parameter space by computing the (log) likelihood 
(usually in the form of ^-statistic) on the uniformly distributed grid points (often referred 
as a template bank). A particular implementation of this method is presented in this paper. 
Another version of the grid based method is described in 18]. There the authors adopt 
the software used for searching continuous GW signals in the LIGO/VIRGO/GEO600 data 
to construct the grid in the parameter space (in four Doppler parameters: sky location, 
frequency of GW signals and its derivative at some fiducial time) and to compute J 7 - statistic. 
Besides that, the authors in 18] attempt to detect several signals at once while we are 
dealing with one (the strongest) signal at the time. They separate the secondary maxima 
from the primary by requiring the coincidence in the parameters recovered from the anal ysis 
of different time delay interferometry (TDI) streams (for more information on TDI, see 19 1 
and references therein). Once the prime maxima are identified, the grid mesh is refined 
by zooming in onto them to improve the parameter estimation. The detected signals are 
removed from the data and the procedure is repeated. 

The last type of search uses Radon transform to identify the frequency of the signal and 
the source's sky location 20] . The basis of this method is that the LISA response function 
can be seen as a Radon transform of binary distribution in those three parameters. 

Galactic binaries were present in all four challenges conducted so far. In this article we 
report our results of the analysis of challenge 3.1 data set. This data set contains approx- 
imately 6 x 10 7 Galactic binaries with simulated instrumental noise. The GW signals had 
measurable frequency evolution at high frequencies. The participants of the challenge have 
returned the parameters of the detected signals: the sky position in ecliptic coordinates, fre- 
quency of GW and its first derivative, inclination of the orbit to the line of sight, polarization 
angle, initial GW phase and the amplitude. 

The paper is organized as follows. In Section 2 we shall present an analytic approximation 
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to the response of the LISA detector to a gravitational wave signal from a binary system. 



We shall not give details of the derivation, these can be found in the original papers 



2l|, 



221 ]. and [23|. In Section 3 we shall present the maximum likelihood method in application 



to detection and estimation of parameters of the GW signal from a binary system imbedded 
in stationary Gaussian noise. In Section 4 we describe our data analysis tools and algorithms 
that we have implemented our computer codes for challenge 3.1 search. We describe the 
search strategy in Section 5 and discuss the results of the search in Section 6. 



II. RESPONSE OF THE LISA DETECTOR TO THE GRAVITATIONAL- WAVE 
SIGNAL FROM A BINARY SYSTEM 

The LISA detector consists of 3 satellites forming a constellation of an approximately 
equilateral triangle. The constellation rotates around the Sun with a period of 1 year trailing 
the Earth by 20 degrees. The triangular constellation is inclined at 60 degrees to the ecliptic 
and rotates itself around its center with a period of 1 year in the direction opposite to the 
rotation around the Sun. The LISA detector in general will produce 3 independent data 
streams. In the long wavelength approximation, when the length of the gravitational wave 
is much longer than the distance between the spacecraft, the number of independent data 
streams degenerates into 2. There are various combinations of the responses of the LISA 



detector (see [21] for details). The data simulated for the mock LISA data challenge are 
the first-generation TDI Michelson combinations that we denote by X, Y, and Z. TDI is a 
software technique (jl]]) by which we remove the dominant frequency noise from the LISA 
instrumental noise, first generation means that in the TDI procedure we assume that the 
distances between the spacecraft are constant, independent of time. In this Section we shall 
summarize approximate analytic formulas for the first TDI generation Michelson responses 
of the LISA detector to a gravitational wave signal from a binary system. The formulas are 



essentially the same as in Appendix C 23( except that they are given using the conventions 



used in the Synthetic LISA numerical software 



24j . The GW response of the first-generation 



TDI Michelson observable X is given by a linear combination of the four time-dependent 
functions X^ k \t). 

4 

X{t)=2uL sm{uL)J2a (k) X {k \t), (1) 



k=l 
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where u is the angular gravitational wave frequency and L is the distance between the 
spacecrafts of the LISA detector. 

Before we give functions X^ k '(t) let us introduce few notations. First we define the 



polarization basis in the solar system barycenter frame (following 25||): 



dk 

u = {sin(/3),cos(A),sin(/3)} ~ — (2) 

r)k 

v = { s in(A),-cos(A),0}~^, (3) 

where (3 and A are, respectively, the latitude and the longitude of the source in ecliptic coor- 
dinates and k = — {cos(/3) cos(A), cos(/3) sin(A), sin(/3)} is direction of the wave propagation. 
Then we introduce the LISA motion used in the production of MLDC data. The position 
of each spacecraft can be split in the position of the guiding center R and position of the 
spacecraft with respect to that center: 

n = R + Lqi, 2 = 1,2,3, (4) 

here we have assumed LISA to be a rigid equilateral triangle: L = L\ = L 2 = L 3 , and the 
qi are as follows 



0% = t-4= jcos(2fit - Xi) ~ 3cos(xi),sin(2f2t - Xi) -3siii(xi), -v / 12cos(fit - Xi)} » ( 5 ) 



2V12 

where Xi = 2(i — l)7r/3. The unit vectors along the arms can be defined via vectors qi: 
^i — 02 ~ 03; an d others are obtained by cyclic permutation of indices: 1 — > 2 — > 3 — > 1. 
Now we are ready to write X^: 



X 0) 



u 2 {t) 
v 2 (t) 



{sinc[(l + kn 2 )x/2] cos[0(t) + (x/2)kq 2 - 3x/2] 



+ sinc[(l - kn 2 )x/2] cos [</>(£) + (x/2)kq 2 - 5ar/2] } 



u 3 (t) 



{sinc[(l + kn 3 )x/2] cos[0(t) + (x/2)kq 3 - 5x/2] 



+ sine [(1 - kn 3 )x/2] cos[0(t) + (x/2)kq 3 - 3x/2] } 



(6) 
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and X^ 3 \ X^ are obtained by replacing cos with sin in and X^ respectively. In the 
above equation we have used: 

Ui = -- [{uhif - (vni) 2 ] (7) 

Vi = (uhi)(vhi), (8) 



where x = loL and sinc(. . .) = The GW responses for Y and Z can be obtained by 

cyclical permutation of the spacecraft indices. 
The phase modulation function <p(t) is given by 

<p(t) =ut+ ^djt 2 + (u + ut) R cos cos(fit + r] - A) , (9) 

where Q = 27r/lyear, rj is the position of the constellation on the orbit around the Sun at 
time t = 0, and R is 1 astronomical unit. The small eccentricity of the Earth orbit (e = 
0.017) can be neglected because it contributes less than one cycle to the phase <f>(t) above 
for gravitational wave frequencies from 0.1 mHz to 12 mHz for which we analyze the data. 
The parameter u is the frequency drift which may occur either due to the gravitational 
radiation reaction or as a result of the tidal interaction between the components of the 
binary system. In the case of a detached binary system evolving only due to the gravitational 



radiation reaction the frequency drift u is approximately given by (see Section IID of 23 1 
for discussion). 

where M. c = vti 2 / (mi + m 2 ) 1 ^ 5 is the chirp mass (mi and m 2 are the individual masses 
of the components of the binary). 

Finally the constant amplitudes take the form 

= Hq cos0o cos2-0 — sin0 o sin 2^, (11) 

= cos 0o sin 2-0 + sin0 o cos 2^, (12) 

= — Hq sin O cos 2ip — cos0 o sin 2^, (13) 

= — Hq sin 0o sin 2-0 + Hq cos0o cos 2ip, (14) 



where 



h+ = ho(l + cos 2 0/2, (15) 
= h cosi. (16) 



The parameters ho, <po, if), and i are constant amplitude, the constant phase of the signal, the 
polarization angle, and the inclination angle respectively. In the case of a detached binary 
system evolving only due to the gravitational radiation reaction the constant amplitude h 
is given by 



ho 



4(OM c ) 5 / 3 run 2/3 



(17) 



c 4 D L 

where Dl is the luminosity distance to the source. 

One can invert the equations (fTTj) for amplitudes to obtain formulas for astrophysical 
parameters h , (f> , if>, and i. We first introduce the quantities 

A = (a«) 2 + (a^) 2 + (a^) 2 + (a^) 2 , (18) 

D = aMa^-aWa®. (19) 

Then the constants h$, h~$ , ho, and <po can be uniquely determined. 



h+ = J(A + VA 2 -4D*)/2, (20) 



/i x = sign(D)J(A - VA 2 - W 2 )/2, (21) 



ho = K + \JK 2 ~ K 2 , (22) 
i = acos(/iQ /h ). (23) 

Finally the constant phase (f>o and the polarization angle if> can be obtained from the following 
equations: 

2 ( a (l)q(3) +Q (2) a (4) ) 

tan20 o - (fl(3))2 + (g(4))2 _ (fl(1))2 _ (a(2))2 , (24) 

2(Q(D a (2) +Q (3) Q (4) ) 

tan 4^ - (a(1))2 + (q(3))2 _ (fl(2))2 _ (q(4))2 . (25) 

The long-wavelength (LW) approximation to the GW responses is obtained by taking the 
leading-order terms of the generic expressions in the limit of uL — > 0: 

X LW (t) ~ 4{uL) 2 {[u 2 (t) - u 3 (t)} [a {1) cos <f){t) + a (3) sin0(t)] (26) 
+[v 2 (t) - v 3 (t)} [a (2) cos0(t) + a (4) sin0(t)]} (27) 

and Y LW , Z LW responses are obtained by cyclical permutation of the indices. 

In Fig. 1 we have compared the power spectra of the signal generated using the analytic 
formulas given above by equation ((TJ) with the power spectrum of noise free response gener- 
ated by Synthetic LISA software for one of the training sets provided in MLDC Synthetic 
LISA software was used to generate the analyzed MLDC data sets. 
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FIG. 1: Comparison of the power spectra of the gravitational- wave signal response using the 
analytic formulas presented in this paper and Synthetic LISA 24]. 



III. MAXIMUM LIKELIHOOD DETECTION AND PARAMETER ESTIMATION 



To detect the signal and estimate its parameters we use the maximum likelihood (ML) 
estimation method which consists of maximizing the likelihood function A with respect to 
the parameters of the signal 26[. Let us first consider the Michelson X combination given 
in the previous section. We assume that the noise in the detector is a stationary Gaussian 
random process. Moreover we assume that over the bandwidth of the signal the spectral 
density of the noise is approximately constant and equal to S a = S(cj ). This condition 
should be well fulfilled for the case of a gravitational wave signal from a white dwarf binary 
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in LISA detector noise. Then we can approximate the log likelihood function by 

logA = 2|[(^X)-i(X 2 )], (28) 

where y x , is the noisy data in the X channel, S Q is one-sided spectral density of the noise, 
T Q is the observation time, and the time-averaging operator (■) is defined by 

(<?>:= ^ f° g{t)dt. (29) 

2 o Jo 

By introducing the variables 

h (k) = 2uL sm{uL)X {k) (30) 
the X combination can be written in a compact form 

X = Y,a^h^ k \ (31) 
k=i 

The ML estimators a^ k ' of the amplitudes are found by maximizing log A with respect to 
parameters cS k \ that is by solving 

The equations (1321) above are equivalent to the following set of linear equations 

4 

^ M (0W a W = iV (0 ) / = !,..., 4, (33) 
k=i 

where 

M (i)(k) = ( h (k) h (i)^ (34) 
= (y*/j(Q). (35) 



Thus the maximum likelihood estimators of the amplitudes are explicitly given by 

a( fc ) = ^(M- 1 ) (0(fc) iV«. (36) 
i=i 

Substituting the above estimators for amplitudes in the log likelihood function logA 
yields the reduced log likelihood function that we denote by T: 

- F = SEE( M " 1 ) (/)(fe)iv(0iv(fc) - ( 3? ) 

1=1 k=l 
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We call the above function the ^-statistic. 

One can show that the following relations hold approximately for the components of the 
matrix M®W 



(h<»h®) = (h^h^) = 0, 
(h<»hW) = (h®h®), 
(h®h®) = (h^h®), 
(h^h®) = (h®h&), 
(h<»hW) = -(h&>h®). 
It is convenient to introduce the following variables 

U = 2(h w h {1) ), 
V = 2(h® /i (2) >, 
Q = 2(hW h^), 
P = 2(/i (1) /i (4) ). 

Let us introduce the complex amplitude parameters 



a« W 3) , 
aW + ia&, 



(38) 
(39) 
(40) 
(41) 
(42) 

(43) 
(44) 
(45) 
(46) 

(47) 
(48) 



where a^ k ' are given by Eqs. (|lip - (|14p . Let us also define the complex modulation functions 
and mS v > by 

{sine [(1 + kn 2 )x/2~\ exp i [(x/2)kq 2 — 3x/2] 

+ sinc[(l — kn 2 )x/2~\ expi[(x/2)kq 2 — 5x/2] } — 
{sine [(1 + +kn 3 )x/2~\ exp i [(x/2)kq 3 — 5x/2] 



m (u) 




u 2 {t) 






y 2 (t)_ 



(49) 



+ sine [(1 - kn 3 )x/2] expi[(x/2)kq 3 -3x/2]}. 
Introducing further complex quantities 

W = Q + iP, 

N (u) = N m + iN (s)^ 

N (v) = N W +iN W 



(50) 
(51) 
(52) 
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where N^ k \ Q and P are given by Eqs. f[3"5"j) . (J15|) . and (T4"6"l) respectively and using the 
approximate relations given by Eqs. ( 1381) . ( 1391) . ( 1401) . (HIT) , and ()42l) ) we can write the ML 
amplitude estimators of the complex amplitudes and the .F-statistic in the following compact 
form 

'aW \ 2 ( V -W* \ ( 
I A \-W U I \ 

|2 . rT \ , t /„m2 



(53) 



T 



V\NV>\ + U\N<?>\ - 2Re (N^ ) *] 



r=2 i* a '- < 54 > 

where A = UV — \ W\ 2 . The integrals and can be expressed as 

N {u) = 2uL sm(uL)(y x {t)m {u \t) expi<f>(t)), (55) 
iV w = 2 sin(wL)(y x (t)m^(t) expi0(t)). (56) 

As we shall see in the next section the above form of the ^-statistic is very suitable for 
a numerical implementation. In a similar manner one can derive the ^-statistic for other 
Michelson variables. 

In order to extract information about the signal from all three independent variables we 
need to derive the J-"-statistic for the whole LISA network. It is then useful to consider the so 
called "optimal" combinations of the responses. These combinations have the property that 
their instrumental noises are uncorrelated (see [27]) and consequently their cross-spectrum 
matrix is diagonal. In this case the log likelihood function for the whole network is the sum 
of log likelihood functions for the individual combinations, For the Michelson variables the 



optimal combination are given by 28] 



A = ^f, (57) 

The noisy data, y A ,y E and y T are obtained as the analogous combination of the Michelson 
observables y x , y Y and y z , where y Y , y z are data in the Y and Z channels respectively. The 
log likelihood function for the network takes the form 

log£ = 2T Q I-*—[(y*A) - l(A 2 }]+ 

{Sa{Uo) 2 , (60) 

' i(y E E)-l(E*)] + -^[(y T T)- 1 -(T>)] 



SeM"" 2 x n S t (u )"» 1 2 
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(62) 



The power spectral densities S A , Se, St are given by 

S A (u) = S E (cu) = 32cos(a;L/2) 2 sin(wL/2) 2 {[6 + 4cos(a;L) + 

2 cos {2uL)]S pm + [2 + cos (loL)}S° p } (61) 
St(w) = 128 cos (loL/2) 2 sin (wL/2) 4 [4 sin (coL/2) 2 S pm + S op ] , 

where S pm and S" 515 are spectral densities of proof-mass noise and optical path noise respec- 
tively. It turns out that the ML estimators of the amplitudes and the ^-statistic can be 
written in the same form as for a single detector. To do this we introduce the following 
noise-weighted average procedure. For any two vectorial quantities p and q, 

p(t) = {p A (t),p E (t),p T (t)) , 
q(t) = (q A (t),q E (t),q T (t)), 

the noise-weighted average operator (-) s is defined as follows, 

(p q) 5 := w A (p A q A ) + w E (p E q E ) + w T (p T q T ) , (63) 
where the weights wj (I = A, E, T) are defined by 

<jT-l 

wj := I = A,E, T, with S- 1 = S A ' + Sz 1 + S T \ (64) 

With the above definitions the log likelihood function of Eq. (J60l) defined for data 
y = (y A y y E \ V T ) an d response R = (A, E, T) can be written in a compact form 

log£ = 2-^[( y R) 5 -i(R 2 > 5 ]. (65) 
Having; now defined by [see Eq. IjBI] )] 

4 



R = ^/)hW (66) 



fc=i 



and A^C*)W, J\fW by [see Eqs. (03), (05])] 

M W) = (hWh») s , (67) 
= (yh«) s . (68) 

it is straightforward to get the maximum likelihood estimators of the complex amplitudes 
and the ^-statistic for the LISA network 



a opt 
a opt 




(69) 
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= 2^- u — (70) 

where A = UV - |W| 2 and W = Q + iP. The quantities V, U, Q, P in Eqs. ([69j and (170]) 
above are defined in the same way as quantities V, U, Q, P in Eqs. (|43l) - fj46|) . but with the 
time-averaging operator (•) replaced everywhere by noise-weighted averaging operator (-) s 
and scalar functions by their vectorial counterparts h^. The two functions Af^ and 
ATM are explicitly given by 

A/» = A/" (1) + iW (3) = 2uL sin(a;L)(y(t)m( u )(t) expi<j)(t)) s , 
ATM = ATW+tAfW = 2uL sm(uL)(y(t)m^(t) expz</>(t)) s , 

where two vector functions and are the relevant combinations of rn^'s and 
defined for X, Y and Z [see Eq. flM}]. 



IV. DATA ANALYSIS ALGORITHMS 
A. Use of the FFT algorithm 

The detection statistic, T [Eq. f )54|) or (170]) ]. involves integrals and of the form 

1= / y(t) m(t;u, /3, X) exp[i(f) mo d(t;u,u, /3, X)]exp[iut] dt (71) 
</o 

where m is one of the two complex modulation functions ( 149]) . while the phase modulation 
(p mod is given by 

0mod(£; /?, X) = -^jjt 2 + uRcos (3 cos(Qt + ?7o — A) (72) 

[see Eq. ([9j)]. In order to evaluate the integral (ITT]) efficiently we would like to use the FFT 
algorithm. The integral ( ITT]) is not a Fourier transform because both the phase modulation 
function mo d and the amplitude modulation function m depend on the angular frequency oj. 
We overcome these problems in the following way. Firstly we introduce a new representation 
of the phase, namely we introduce two new parameters 

A = uRcos (3 cos(A — 770) , 

B = uRcos (3 sm(X — r] ). (73) 
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In this new parametrization the phase modulation takes the form 

0modO; w, A, B) = \jot 2 + A cos(fit) + B sin(fit), (74) 
and it is independent of the angular frequency parameter u. A similar parametrization 



has been used in the search of the resonant bar NAUT 



gravitational waves from spinning neutron stars (see 29j). (As we shall see in the fol- 



LUS detector detector data for 



lowing section the above parametrization by the coordinates {u, u, A, B} will also prove 
useful in the construction of our grid.) Secondly we assume that the bandwidth [w! w 2 ] 
of the data is small so that the amplitude modulation function m varies little over the 
interval of angular frequency from lower edge of the band u)\ to the upper edge of the 
band Co>2- In order to satisfy the above approximation, in our search we have divided the 
data into narrow bands of Af = {0J2 — coi)/(2ir) = 0.1 mHz by passing them through 
narrowband filters (see Section |V] for details) and each narrow-banded data were an- 
alyzed separately. Then we can approximate the modulation function m(t; u, (3, A) by 
m(t;uj mid ,A,B) = m(t; ui mid , (3 mid (A, B, ui mid ), \ mid (A, B, u mid )) where uj mid = (u 2 - Wi)/2 
and where (3 mid (A, B, u mid ) and X mid (A, B , oj m i d ) are obtained by inverting Eqs. (173|) for 
lo = uj mid . They are explicitly given by 



(J j\2 _|_ £>2 \ 
5— > ( 75 ) 

A mid (A, J B,a; mid ) = i] + arctan . (76) 

Note that there are two values of the ecliptic latitude (3 for each pair of the values A and B. 
Consequently the integral ( |7T1) can be approximated by 

X~ / y(t)m(t;u mid ,A,B)exp[i(p mod (t;uj,A,B)]exp[iujt]dt, (77) 
^0 

For discrete data y(t) the above integral can be converted to a discrete Fourier transform 
which can be calculated by the FFT algorithm. 



B. Metric on the intrinsic parameters space 

In our method the detection of weak, quasi-monochromatic GW signals relies on an 
efficient placement of the templates in the bank. It should minimize the number of templates 
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for a certain accepted loss of a signal-to-noise ratio. Here we follow the geometric approach 



initialized in 
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3 11 ] and introduce a metric on the intrinsic parameters space to measure the 



mismatch between the true signal and the template. In order to construct the metric and a 
grid on the intrinsic parameter space over which we calculate the ./-"-statistic we introduce 
an approximation to the signal that we call the linear model [32]: 

1 



s(t) = A cos (cot + -cut + A cos Vtt + B sin Vtt + O ) 



(78) 



where A is a constant amplitude and 0o is a constant phase and where A and B are parame- 
ters given by Eqs. (173]) . The phase modulation of the linear model is exactly the same as that 
of the exact model whereas the amplitude is constant. This is a reasonable approximation 
because the amplitude modulation functions vary very slowly, they are periodic functions of 
one year. 

The metric on intrinsic parameter space is defined by the reduced Fisher matrix which 
is obtained from the full Fisher matrix of the linear model by projecting on the intrinsic 
parameters space and normalizing it. The reduced Fisher matrix V determines the loss 
of signal-to-noise ratio when parameters of the signal, 9 = (cu, Cj, A, B), differ from the 



parameters of the template by A9 = (Acu, Acu, AA, AB) 33j, |34J ]: 



T e iA9,A9) 



p|(q) - p}m 



+ 0(\A9\ 



(79) 



where pg{0) is the optimal signal-to-noise ratio and p^(A9) is the signal-to-noise ratio for 
the mismatch A9. The quantity r 2 is usually called the mismatch and it is denoted by m 
in 341 ] . For the calculations of the reduced Fisher matrix we refer the reader to Appendix 
lAl We show there that the approximation of the linear model leads to a particularly simple 
forms of the reduced Fisher matrices on 3-dimensional intrinsic parameter space {u, A, B} 

0i() 



2irn 



\-2^ 



(80) 
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and 4- dimensional intrinsic parameter space {u, u, A, B} (lAllj) 
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We will use the reduced Fisher matrices ( 150]) and f lHTj) to build the grid in such a way 
that the distance defined by 1^ or T from any point of the parameter space to the nearest 
node of the grid is not larger than some fixed value r. We also see that not only the metrics 
T 3 and T are flat but in the coordinates {u, A, B} and {u, u, A, B} their coefficients are 
constant, independent of the values of the parameters. 



C. Construction of the grid in the parameter space 

1. Covering problem on lattices 



The problem of constructing a grid in a d- dimensional parameter space is equivalent 
to the problem of covering (^-dimensional space with equal overlapping spheres of a given 
radius. The optimal covering would have minimal possible thickness or density of covering 
defined as the average number of spheres that contain a point of the space (see 35| and 
below). When the metric is flat, as in the linear model, centers of the spheres can lie on a 
c?-dimensional lattice. In this case one can take advantage of the theory of lattice coverings. 
In the rest of the chapter we briefly sketch the basic definitions from the theory of lattices 
that will be used in the construction of the grid. 

In general for any discrete set of points S = s*2, . . .} in MJ 1 the covering radius R of 
S is defined as the least upper bound for any point of M. n to the closest point sf. 

R(S) = sup inf \x — s\. 



ses 



Then spheres of equal radius r centered at the points Si will cover R n only if r > R. 

A lattice A is a discrete subset of M. n . Any lattice has a basis b = j&i, b n \ of linearly 
independent vectors on K n such that the lattice is the set of all linear combinations of fej's 
with integer coefficients: 



A 



^ c$i : Cj G Z, 



1,2, 



n 



(82) 
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A lattice basis is not unique, in dimensions d > 1 there are infinitely many of them, but all 
the bases have the same number of elements called the dimension of the lattice. To specify 
a basis b of a lattice we will use the notation A(6). 

A lattice Ai is equivalent to a lattice A 2 if Ai can be transformed into A 2 by a rotation 
reflection and change of scale. 

The parallelotope consisting of points c\b\ + . . . + c n bi with < q < 1 is a fundamental 
parallelotope and is an example of an elementary cell, that is the building block containing 
one lattice point which tiles the whole K n by translations of lattice vectors. There are 
infinitely many elementary cells but the volume of each elementary cell is unique for a given 
lattice A. 

The Voronoi cell around any point v of A is the set of vectors x of M. n which are closer 
to v than to any other lattice vector: 

V(v) — {x : \x — w\ > \x — v\ for all w G A} . (83) 

All Voronoi cells of a given lattice are congruent convex polytopes and are another examples 
of elementary cells sometimes referred to as Wigner-Seitz cells or Brillouin zones. 

For the lattice A having Voronoi cells congruent to polytope V(v), where v is any of the 
lattice points, the covering radius R(A) is the circumradius of V(v) i.e. the largest distance 
between v and the vertices of V(v). 

The thickness of the lattice covering is given by 

volume of <i- dimensional sphere of radius R(A) 
volume of the elementary cell of A 

The covering problem asks to find a lattice with the lowest thickness. The thinnest lattice 
coverings are known in dimensions up to 5. They are given by the so called Voronoi's 
principal lattices of the first type and are denoted by A* d . A* 2 is equivalent to the hexagonal 
lattice and is proved to be thinnest covering of the plane, A\ is equivalent to the body- 
centered- cubic (bcc) lattice. For the results of the best known coverings in higher dimensions 



we refer readers to 



35]. 



2. Covering problem with constraints 

In our search scheme the calculation of the ^-statistic involves two Fourier transforms 
that can be computed efficiently using the fast Fourier transform algorithm. For this reason 
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we want the nodes of the grid to coincide with Fourier frequencies: Au, 2Au, 3Au, ... for 
some fixed frequency resolution Au. This imposes a condition that one of the lattice basis 
vectors has a fixed length 



and forbids an immediate use of the general results of the theory of lattice coverings. Instead 
one can formulate the covering problem with constraint: to find the thinnest lattice covering 
of the d- dimensional space with spheres of radius r and one of the basis vectors of the lattice 
having fixed length As far as we know the general solution to the problem is not known. 
We present a construction of a nearly optimal lattice that satisfies the constraint with a 
good accuracy. 

Let a vector v define the frequency resolution. We search for a lattice A(w') of covering 
radius R(A(w')) = r with lattice basis w' satisfying the constraints that can be expressed as 
w' = {vo, w[, . . . , We find the thinnest constrained lattice starting with an optimal 

unconstrained lattice in ^-dimensions. The idea is to shrink the optimal lattice as little 
as possible such that one of the basis vectors of the resulting lattice coincides with the 
constraint vector v . We notice that the orientation of the constraint vector v has no effect 
on the optimal constrained lattice, it is only the length of vq and assumed value of covering 
radius that matters (and more precisely: only their ratio because the overall scale can be 
taken arbitrarily). 

For a given lattice A there always exists the lattice vector I such that ||v | — |/|| takes 
minimum value that we denote by 1(A). We define Algorithm 1: 

Listing 1. "Minimal" deformation of a lattice. One of the nodes 
of the final lattice coincides with the resolution vector 
Input: Lattice A; vector v . 
Output: Lattice A'; f(A') = v 
Al. Find 1(A). 

A2. Contract A along 1(A) to obtain A c with |Z(A C )| = |ub|- 
A3. Rotate A c to obtain a lattice A rc with l(A rc ) = v . 
A4. Return A' = A rc . 
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For optimal initial lattices Algorithm 1 defines the following function f$ : R — > R: for 
x = R(A), fv (x) = R(A'). For a given r G R we denote by the value of x for which 
the function \f$ (x) — r\ reaches its minimum. The optimal constrained lattice is obtained 
by application of Algorithm 1 to an optimal (unconstrained) lattice A with covering radius 
R(A) = r t . 

In dimensions d = 2, 3, 4, 5 as the initial lattices one takes A* d lattices but the procedure 
can be generalized to any number of dimensions by taking as the input the best known lat- 
tice covering in a given dimension 35] . Fig. [2] illustrates the procedure for two-dimensional 
A* 2 hexagonal lattice. It is seen there that contraction of the lattice can change the initial 
Voronoi cell and covering radius. However when the vector 1(A) initially lies on the "con- 




v =l(A rc ) 




FIG. 2: Illustration of Algorithm 1 in two dimensions. Initial optimal lattice (top left) with the 
basis {wi,W2}, resolution vector vq defining constraint surface and vector I is contracted along the 
vector I (top right) and rotated (bottom left) such that for the new lattice I = vq. Final suboptimal 
lattice has basis {vqjW^} (bottom right). 



straint surface" depicted by the dashed circle in the Fig. [2] the procedure acts trivially (no 
contraction only rotation) leaving the final lattice optimal. One often encounters trivial 
procedures when the vector 1(A) is large as compared to the Voronoi cell and moreover in 
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these cases contractions are small. On the other hand the last trivial procedure occurs when 
the resolution vector and the shortest lattice vector have equal lengths. 

One can determine the values of covering radii for which the final lattices are opti- 
mal. As an example we find the values of four largest covering radii having this prop- 
erty. To do this we consider the sequence of the shortest vectors of the lattice. In three 
and four dimensions the first nonzero shortest vectors of A 3 and A\ lattices have lengths 
2y^3/5R(A* 3 ), A/V5R(A* 3 ), 4t/2/5R(A* 3 ), 2 v /1/1LR(A*), ... and y/2R(Al), y/3R(AQ, 
\/ER(Al), \nR(Al), . . . respectively. The length of the resolution vector which in dimension- 
less units has the from v = (2tc, 0, . . .) is equal to ^/f(t/ , Vq) = tc/V3 in both dimensions. 
This gives the following squares of the largest covering radii for optimal constrained lattices: 
R(A* 3 ) 2 : 5/36tt 2 ps 1.3708, 5/48tt 2 ps 1.0281, 5/96tt 2 ps 0.514, 5/132tt 2 ps 0.3739, ... and 
R(Al) 2 : tt 2 /6 ~ 1.6450, vr 2 /9 « 1.0966, vr 2 /15 ps 0.6580, tt 2 /21 ps 0.4700 [points a, b, c, d on 
Fig-©]. For larger values of R thickness of final lattices grows monotonically with covering 
radius. 

These features are seen in the Fig. [3] which shows the results of application of the 
construction to the model ( lAlj) in 3 and 4 dimensions for different covering radii. The 
resolution vector is Vq = (2ir, 0, . . .) and the observation time is 2 years. As for a fixed 
volume of the parameter space the number of points needed to cover the space with balls of 
a given radius (i.e. for allowed loss of signal-to-noise ratio) is proportional to the thickness, 
the diagram demonstrates that the excess of points due to constraints is minor for the wide 
range of radii which makes the search strategy based on FFT effective. 



D. Number of false alarms 



In order to estimate the number of false alarms expected in our search we use a general 
approach consisting of dividing the parameter space into elementary cells defined by the 



autocorrelation function of the J-"-statistic 



331 ] and |36| Chapter 6.1.3). We use the Taylor 



expansion of the autocorrelation function around the true values of the parameters 9k and 
moreover we use an approximate response of the detector given by the linear model (Eq. 
( 1A14jl or (1A15I) ). In this case the hypervolume of the elementary cell is given by the volume 
V c of the hyperellipsoid defined as 



f H A0 fc A0, < 1/2. 



(85) 
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FIG. 3: Thickness of the lattices in 3 and 4 dimensions compared to the optimal thickness of A% ~ 
1.4635 and A\ ~ 1.7655. The upper diagram shows also the optimal lattice with covering radius 
R(A* 3 ) 2 = 5/48vr 2 . The Voronoi cell, basis vectors, resolution vector and part of the constraint 
surface are depicted; for this specific value of R the head of the resolution vector lies on the 
constraint surface. 



Thus V c is given by 



m/2 



(86) 



T(m/2 + 1) y/detf 

where m is the dimension of the parameter space and V denotes the Gamma function. The 
determinants of the reduced Fisher matrix in 3 and 4 dimensional cases read 

T 2 n 2 n 2 - 6 



det T, 



48 7r 2 n 2 



(87) 
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where it is assumed that the observation time T is an integer multiple n of years. As for 
the linear model the reduced Fisher matrix has components independent of the values of the 
parameters and the number of cells is simply given by 

iV c = (89) 

where V is the hypervolume of the parameter space. We assume that we search a narrow 
band of bandwidth Au with upper frequency w max . Then the volume of the parameter space 
V3 when the frequency drift is not included in the search is given by 

V 3 = 2Aunu 2 max R 2 (90) 

whereas the hypervolume V4 when it is included is given by 

V A = V 3 Au, (91) 

where Au is the range of the frequency drift parameter u. The factor of 2 in Eq. ( 190]) is 
because we search the space spanned by parameters u, u, A, and B twice - both for positive 
and negative values of the ecliptic latitude (3. 

The expected number of the false alarms Np is given by 

N F = N C P F {7 ), (92) 

where Pp is the probability of false alarm. In the case of linear model Pp is the x 2 probability 
distribution with two degrees of freedom i.e. 

Pf(T ) = exp(-jr o ) (93) 



E. Computation of the ^-statistic 

The ^-statistic is computed approximately taking advantage of the speed of the FFT 
algorithm as described in Section IIV Al The J 7 is calculated on the grid constructed in 
Section IIV CI For each parameter pair (A, B) the ^-statistic is computed for both positive 
and negative value of the ecliptic latitude /3 mi d (see Eq. f J75|) ). Approximate calculation of the 
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^-statistic on the grid described above is called the coarse search. Using the coarse search 
we identify signals for which the ^-statistic crosses a certain threshold. The coarse search is 
then followed by the second step which we call the fine search. Fine search consists of a search 
for the maximum of the ^-statistic around the parameters identified by the coarse search. 
To find the maximum we use the Nelder-Mead maximization algorithm 37]]. In the fine step 
we use accurate expressions for the ^-statistic [Eq. ( 154]) or ( !70|) ]. without approximations 
described in Section IIV Al As initial values for the Nelder-Mead algorithm we use the 
parameters obtained in the coarse search. The values of the parameters corresponding to 
the maximum of ^-statistic are our final estimates of the parameters of the signal. 



V. SEARCH STRATEGY 



In our entry for challenge 3.1 we have used the following procedure to extract GW signals 
from white dwarf /white dwarf binaries in the mock LISA data. We search the band from 
frequency / = 0.1 mUz to frequency 12 mHz where / = u/2ir. We do not go to higher 
frequencies because with our search strategy it would involve much more computing time 
than we could afford. Also above 12 mHz the expected number of white dwarf binaries in 
our Galaxy is very small. We first divide the data into bands of 0.1 mHz each. To obtain 
narrowband data in the frequency band [f\ f?\ we first pass the data through 3rd order 
Butterworth filter with passband of [f\ — e f'2 + e] where we choose the edge parameter e 
equal to 0.005 mHz. Then we shift the data to DC by frequency fi — e and we pass it 
again through the 3rd order Butterworth filter with passband of [f% — fi + 2e] . After each 
Butterworth filter we downsample the data. This reduces the number of data points by a 
factor of around 300. Each Butterworth filter is applied twice: forward and backward in 
time. In this way there is no phase shift of the narrowbanded data with respect to the 
original one. In the search of each narrow band for signals we neglect the edges of the band 
of e = 0.005 mHz and therefore we only search the band [fx f?\. We have included in our 
search the frequency drift parameter u. Analysis of the accuracy of estimation of u have 
shown (see Section I VI I below) that it is useful to include the Cj parameter in the search for 
frequencies above 3mHz. We have selected the range of the u parameter using a fit from 
the values of the u parameter in the key of the challenge 3.1 training data set. In each band 
we search for the signals calculating the .F-statistic over the constrained grid constructed 
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in Section IIVCI This enables application of the FFT algorithm in F-statistic computation. 
We select the strongest signal and to this signal we apply the fine search described in Section 
HVEI to estimate its parameters. We use the Nelder-Mead algorithm with the initial values 
provided by the parameters of the template with the largest value of the F-statistic over 
the grid. We reconstruct the signal in the time domain and remove it from the data. We 
then search for the next strongest signal and so on until the signal-to-noise ratio (SNR) of 
the detected signal estimated as 

SNR = y/2(F - 2) (94) 

falls below a certain threshold. We have chosen the threshold for the ./-"-statistic equal to 
18. This corresponds to SNR threshold of around 5.7 (see Eq. (I94j) ). The number of signals 
increases as the signal-to-noise ratio decreases. At a sufficiently low signal-to-noise ratio the 
signals are so close to each other in the parameter space that they interfere with each other. 
The .F-statistic (Eq. l7D|l that we use in our search was derived under the assumption that 
there was only one signal present in the data so it can be used to detect multiple signals 
when they are sufficiently separated in the parameter space so that the /-"-statistic for the 
multiple signals is a sum of the /-"-statistics for individual signals. Thus when there are 
many signals present we neglect the interference between the signals. Our choice of 18 of 
the threshold for the /-"-statistic was a convenient choice that, as we shall see in the following 
section, led to a detection and accurate estimation of over 10 4 signals. 

In the Figs. H] and [5] we have presented application of the above strategy to the challenge 
3.1 data set in the bandwidth from 5 mHz to 5.1 mHz. 

In this bandwidth we have identified 132 signals altogether out of 168 present. The 
accuracy of estimation of most of the signals is one sigma where sigma is calculated from 
the Fisher matrix. For several signals the error was very large indicating that either noise 
mimicked the signal or the residuals of removed signals remained significant. This could 
happen as a result of interference of the signals. 

VI. MLDC RESULTS 

Here we present results of our entry for challenge 3.1 [7|. The challenge 3.1 consisted 
of a two-year data set with 15 s sampling time with signals from around 6 x 10 7 binaries. 
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FIG. 4: Estimation of the signals from white dwarf binaries in the challenge 3.1 data set for the 
band from 5 mHz to 5.1 mHz. Black color denotes the original data and the light blue color is the 
data after signals removed. 

Our aim was to detect as many as possible from the 40628 brightest binaries present in the 
data. We have performed a self-evaluation of our search by the following procedure. In our 
procedure we have used the correlation C(si, S2) between the two signals s\ and S2 denned 
as 



C(si,s 2 ) 



(95) 



y/(si,si) s y/(s2,s 2 ) s 

The first step of our self-evaluation consisted of selection of the true detected signals from 
the set of all submitted data. For a submitted signal s with parameters 8 S we considered set 
B s of all key signals within the frequency bins about f s , i.e. signals satisfying 

\f s - / fc | < l/T obs « 1.6 x 10" 8 Hz, 

where parameters 9k of the key signals were taken from the set of 40628 bright Galactic 
binaries Q. Next, each signal s was paired with the key signal k s from B s that maximizes 
the correlation: 

s 1— > k s : k s = argmaxC (s, k) ; (96) 

fc£B 3 

we interpreted 9 S as the parameters estimates of the key signal k s . In the case of multiple 
detection of a key signal fco, that is when k Bl — k S2 — . . . — ko for different Si, s%, . . ., as the 
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FIG. 5: The zoom of the Fig. [J] around the strongest signal identified, 
main signal Sk we singled out the one that maximizes the correlation with ko, 

s ko = arg max C (s, k ) , (97) 

s: k s =k 

and we rejected the remaining secondary signals. 
A. Original results 

Our original entry for challenge 3.1 contained a bug in our data reading procedure. Our 
analysis for frequencies above 3mHz was performed on the challenge data set shifted in 
time by 1440s. In the original search we have detected 14838 signals altogether and we have 
estimated their parameters. Results of detection are displayed in Fig. [6]which shows detected 
(main) signals and the correlations with signals from the key. We see that the histogram 
of the correlations shows an excess of anticorrelations due to the bug in our data reading 
procedure. Effects of this systematic error can be seen in the Fig. [7] displaying correlations 
as the function of the frequency. Peculiar oscillations of the correlation with frequency are 
easily explained by noticing that the phase difference between nearly monochromatic signal 
and the same signal shifted in time by t Q is equal to 2Trft Q and gives periodically changing 
overlap between the two functions which oscillates with respect to / with "frequency" t 

27 



2000 
1500 
1000 
500 




main signals 
MLDC3 



0.000 0.002 0.004 0.006 0.008 0.010 0.012 
f[mHz] 



1000 
800 
600 
400 
200 





















1 














-1.0 



-0.5 0.0 

correlations 



0.5 



1.0 



FIG. 6: Detection and correlations for the blind challenge 3.1 data set in the original run. The left 
panel shows number of the signals detected by our search and selected by the procedure described 
in the text as a function of the frequency. The right panel displays the histogram of the correlation 
functions (Eq. (]95p ) between our estimated signals and the ones form the key. 
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FIG. 7: Correlations vs. frequency. The periodic variation of the correlation function was due a 
time shift in reading of the challenge data set. 

or by explicitly computing correlation of two monochromatic waves cos [2irf(t — to)] and 
cos(27r/). This explains the excess of signals with anticorrelations in Fig. [6] 
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B. Corrected results 



We have corrected the bug in the reading procedure and have made a second run of the 
challenge 3.1 data set. Here we present results of this second run. In our search we have 
detected 17051 signals altogether and we have estimated their parameters. 

In the band from O.lmHz to 3mHz we have made two runs - one with and the other 
without the Co parameter included. We have used these two runs to investigate at what 
frequency it is useful to start estimating the frequency drift of the GW signal. In the 
search including the Co parameter, we have found that we can start estimating it only above 
frequency of 0.5 mHz because below this frequency the cell of our 4-dimensional grid is 
bigger than our parameter space. Including the frequency derivative from 0.5mHz we have 
identified 16785 signals and including Co starting only from 3mHz we have found 17051 
signals. We have compared the absolute values of the errors in Co when we include it in the 
search and when we do not include it. When we do not include it we take as the error the 
absolute value of the Co parameter from the challenge key signal. We calculate the mean 
values of the these absolute errors in each band. We find that slightly below the frequency 
of 3mHz the mean error for the case when we include Co parameter becomes less than when 
we do not include it. Thus we find that our initial choice of threshold frequency equal to 
3 mHz to include the Co is a reasonable one. Therefore, as our final result, we consider the 
signal found using the search that turns on the frequency derivative at 3mHz frequency. We 
have also estimated the expected number of false alarms using the formula ( 1921) in Section 
IIVDI We find that for low frequencies where we detect most of our signals the number of 
false alarms is negligible. The number of false alarms increases quadratically with frequency 
and linearly with the range of Co parameter. The expected number of false alarms exceeds 
one only at the frequency / = 8 mHz. 

In Fig. [8] we present the number of signals detected against all the challenge 3.1 signals 
and the correlations of the estimated signals with the key signals in the second run. We see 
that there is still an excess of correlations with correlation parameter around zero. We have 
investigated the number of correlations as a function of signal-to-noise ratio (see Fig. [S]). We 
find that the excess of low correlations originates from frequencies below 3 mHz. Moreover 
for SNR > 10 the number of small correlations considerably decreases. The number of 
selected signals after we discarded signals below SNR = 10 and for frequency less than 3 
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FIG. 8: Detection and correlations for the blind challenge 3.1 data set in the second run. 
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FIG. 9: Number of correlations of estimated signals with key signals as a function of the signal-to- 
noise ratio for signals with frequency below 3 mHz. 

mHz is 12805. We quote this number as the number of signals which we detect and accurately 
estimate the parameters. We need to refine our data analysis methods in order to extract 
reliable estimates of the parameters for the signals of SNR less than 10 and frequency less 
than 3 mHz. 

The excess of zero correlation signals arises either because for some low signal-to-noise 
ratio the parameter estimation was not accurate or/and because at low signal-to-noise ratio 
there are much more signals that interfere causing biases in the parameter estimators. Fig. 
[10] demonstrates the results of the parameter estimations for our search of challenge 3.1 data 
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FIG. 10: Errors of the parameters of the signals detected and verified in our search of the challenge 
3.1 data set. 
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FIG. 11: Errors of the estimation of parameters as in Fig. [10] but normalized by the variances 
obtained from the inverse of the Fisher information matrix. 



set. Errors are defined as differences between the key and the recovered signal parameters, 
A8 := Okey — Orec- Histograms in Fig. [TTIshow the parameter estimation errors divided by the 
standard deviations <Jg key obtained from the Fisher matrix for the key signals. In Fig. [12] 
we have presented the power spectrum of the challenge 3.1 data against the power spectrum 
of the data after the detected signals were removed. We plot two power spectra: one when 
all the signals identified are removed and the other one when only the ones selected by our 
procedure are removed form the challenge data set. We notice two effects. One is that 
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le-46 



0.0001 0.001 0.01 

Frequency (Hz) 

FIG. 12: Amplitude square of the Fourier transformed data. Upper signal is original challenge 
data (A-channel, includes full Galaxy and instrumental noise), middle is the data with selected 
removed signals, and lower is the data with all identified signals removed. 

periodically our identification procedure gets worse. This because we were not estimating 
parameters of the signal very well at the edges of the narrow bands. We have found that this 
is because our narrowband filtering procedure was not perfect. The other effect is that when 
we remove all identified signals (lower signal on the Fig. [T2"|) we are doing a much better job 
then when we remove only signals identified as true signals by our procedure (middle signal 
on the Fig. [T2|) . Thus we fit quite a large number of signals well but with wrong parameters. 
In Fig. [13] we have compared smoothed spectrum of the challenge 3.1 data set with that 
of data with identified signals removed and we have compared them with spectrum of the 
LISA detector instrumental noise. From this figure we conclude that above frequency of 6 
mHz we resolve all the white-dwarf binary systems well. 
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FIG. 13: Power spectral density. Upper signal is the original challenge data, middle is reduced 
data (after removing all found signals), and lower is the instrumental noise. 



VII. CONCLUSION 



From the analysis of our challenge 3.1 results we see that in order to increase substantially 
the number of signals with good parameter estimation we need to assume in deriving our 
filters that there are more than one signal present in the data. However even with the current 
procedure that estimates signals one by one we can still make the following improvements. 

1. Improve the splitting of the time series into narrow bands. 

2. Lower the threshold for detection. 

3. For high frequencies bands where signals are well separated identify all the signals at 
one scan of the parameter space instead extracting only one and scanning the whole 
parameter space again to go to the next. 
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Appendix A: Linear model 



In this appendix we consider the linear model, an approximation of the full gravitational 
wave signal introduced in order to construct a metric and a grid on the intrinsic parameter 
space. It has the form: 

1 



s(t) = A cos (cot + + A cos Qt + B sin Qt + 



(Al) 



with constant amplitude A , constant phase 0o an d two parameters A, B (related to u, (3 
and A). 

Following the general case ([I]) we rewrite the signal flAlj) in the form 



s(t) = A 1 X 1 (t) + A 2 X 2 (t), 



(A2) 



where 

Xi = cos (out + — out 2 + A cos Qt + B sin Qt) , 

X 2 = sin (cut + -Cut 2 + A cos Qt + B sin Qt) 

and where we have introduced the extrinsic amplitudes A\ = /i o cos0o an d A 2 = —h sm(j) Q . 
The Fisher matrix, 

(A3) 



2T D 



5"(w ) 

for the linear model (1A2I) with respect to the parameters (A\, A 2 , to, ou, A, B) takes the form 
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(A4) 



where 



Gx 



G 2 



G, 



1 









1 



A 1 2 +A 2 2 
T0A2 T 2 A 2 



2{A X 2 +A 2 2 ) 
_ T0A1 

To 3 



&{A 1 2 +A 2 2 ) 

To 2 A! 
6(Ai 2 +A 2 2 ) 




Jo io Q 



To_ 2 

3 

T 3 



4r 



_1 

Q 



rp 4 

20 

Ik 
2ft 



J_ 

n 2 

1 
2 



20 



(A5) 
(A6) 



(A7) 
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and where the optimal signal-to-noise ratio, 



P 



2T Q 
S(u ) 



(s 2 ), 



(AS) 



is given by p 2 = T °^, + f 2 ^ and the observation time T is assumed to be an integer n of 
years. In the derivation of Eqs. ( 1A5I) - ( 1A8I) we have used the following approximations 



(X, 2 ) ~ (X 2 2 ) 



(X x X 2 ) ~ 



(A9) 



corresponding to approximations fl38|) - fj46|) for the case of the full signal. 

The reduced Fisher matrix obtained from the full Fisher matrix T by projecting V on the 
intrinsic parameters space and normalizing it is explicitly given by (|33j) 

~ 1 



(G3 — G2 T Gi 1 G2) 



(A10) 



The coefficients of the reduced Fisher matrix T of the linear model (IA2I) in the dimen- 
sionless units with T = 1 are given by 
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For a network of detectors s = (sa, se, St) with uncorrelated noises the Fisher matrix 
and the optimal signal-to-noise ratio can be written in terms of the noise-weighted averaging 
operator and vectorial response [see Sect. flHIl) ]: 



opt ij 



Popt 



2T 



■(disdjs) 



3 a ISi 



2T, 



S(u. 



(A12) 



(A13) 



In the case of network of LISA detectors the optimal responses A, E, and T can be approx- 
imated by the linear model of the form 

1 



si(t) = cih Q cos (cut H — cot 2 + A cos Vlt + B sin Vtt + (fi + di), 

2 



I = A,E,T, (A14) 



where constants cj and dj have been introduced in order to take into account different am- 
plitude and phase modulations for each observable. The reduced Fisher matrix for network 
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turns out to be exactly the same as the reduced matrix for a single response. This is the spe 



34J, Ch. IIIC) that grid resolution 



cial case of the general property discovered by R. Prix 
in the parameter space is independent of the number of detectors. 

For low frequencies the frequency derivative Cj is small and there is no need to include 
this parameter in the search. Then the linear model simplifies to 



S3i{t) = ciho cos (cot + A cos fit + B sin fit + <p + dj) , I = A, E, T 



(A15) 



and the corresponding reduced Fisher matrix T3 reads 



( J- x - 

12 u 2im 

o I 



(A16) 
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